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The  effect  of  plasticity  on  the  rate  of  growth  of  fatigue  cracks  is 
significant  for  a  wide  range  of  problems  associated  with  the  damage  tolerance 
assessment  of  aerospace  structures.  The  range  of  problems  includes  crack 
growth  from  cold  worked  fastener  holes,  crack  growth  through  plasticity  due  to 
local  notch  stresses,  crack  driving  force  for  thermal  gradient  fields  and 
welding  residual  strain  fields,  small  flaw  growth  in  high  nominal  stress 
fields,  and  numerous  related  problems.  These  problems  have,  of  course,  been 
analyzed  using  a  variety  of  approximate  analytical  or  numerical  procedures. 
However,  as  will  be  summarized  within  this  report,  many  of  these  earlier 
modeling  approaches  nave  involved  errors  which  may  significantly  affect  the 
predicted  fatigue  crack  growth  life  of  the  structure.  The  current  research 
nas  resulted  in  some  new  and  critical  insights  into  this  class  of  problems, 
while  providing  a  basis  for  improved  modeling  of  these  problems. 

The  current  research  makes  use  of  the  boundary  integral  equation  (BIE) 
method,  as  modified  to  account  exactly  for  the  elastic  crack  problem.  The 
usual  BIE  formulation  for  elastic  problems  reduces  the  numerical  problem  to 
one  of  modeling  the  boundary  data,  while  preserving  the  complete  interior 
solution  of  the  field  equations.  In  the  elastic  fracture  mechanics  problem, 
the  Green’s  function  approach  is  used  wherein  the  BIE  is  modified  to  account 
for  the  presence  of  a  stress  free  crack  at  an  arbitrary  location  within  the 
structure.  The  use  of  the  Green's  function  for  the  crack  eliminates  the  need 
to  model  the  boundary  of  the  crack,  and  provides  a  complete  mathematical 
description  of  the  elastic  strain  field  within  the  body,  due  to  the  crack. 

This  clearly  contrasts  with  the  finite  element  method  which  requires  that  the 
crack  surface  and  the  interior  strains  be  modeled  with  some  set  of 
interpolation  functions. 

The  BIE  method  has  been  successfully  modified  to  account  for 
elastoplastic  response  by  a  number  of  investigators.  However,  extension  of 
the  fracture  mechanics  model  with  the  Green's  function  approach  has  not  been 
previously  demonstrated.  In  order  to  account  for  elastoplastic  response  with 
the  BIE  method  one  must  numerically  model  the  interior  plastic  strain  field. 

In  all  other  ways  the  elastoplasticity  solution  uses  the  standard  elastic  BIE 
formulation.  The  current  work  reports  on  the  successful  extension  of  the 
special  Green's  function  formulation  for  the  fracture  mechanics  problem  to  the 
elastoplasticity  formulation.  Not  only  has  the  work  resulted  in  accurate 
models  of  crack  tip  plasticity  for  a  reference  problem,  but  it  has  shown  some 
important  new  analytical  and  numerical  results  for  cracks  growing  in  plastic 
strain  fields. 

The  second  year  of  the  contract  effort  will  focus  on  the  crack  growth 
problem.  That  is,  the  effect  of  crack  tip  plasticity  on  the  subsequent  crack 
growth  rate  will  be  studied.  The  effects  of  crack  tip  overloads  on  retardation 
or  acceleration  through  closure  and  residual  stress  effects  will  be  studied. 

In  addition,  the  elastoplastic  BIE  formulation  will  be  more  fully  exploited 
for  problems  of  crack  growth  in  residual  strain  and  thermal  strain  fields.  The 
purpose  of  all  of  these  studies  will  be  to  identify  modeling  shortcomings  in 
current  practice  for  these  problems  as  well  as  to  provide  some  new  results  for 
the  small  flaw  problem. 
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RESEARCH  OBJECTIVES 

Generally  speaking,  advanced  aerospace  structures  have  been  designed  for 
damage  tolerance  considerations  using  elastic  fracture  mechanics  models. 
Problems  associated  with  residual  plastic  strains  at  notches,  cold  worked 
fastener  holes,  weld  residual  strains,  and  thermal  gradient  loading  have  been 
modeled  using  elastic  superposition  methods  along  with  elastic  fracture 
mechanics  models.  Crack  tip  plasticity  is  involved  in  all  fatigue  crack 
growth  problems.  Crack  tip  plasticity  dominates  the  problem  of  predicting 
crack  growth  under  spectrum  loading  conditions  where  acceleration  and 
retardation  effects  are  important.  Finally,  the  small  flaw  problem,  wherein 
crack  growth  rate  is  apparently  accelerated  relative  to  the  large  flaw 
problem,  cannot  be  currently  explained  by  elastic  fracture  mechanics 
considerations . 

The  development  of  improved  models  for  the  crack  growth  problem  for  the 
full  range  of  these  proolems  is  crucial  to  the  damage  tolerance  assessment 
needs  for  advanced  aerospeace  structures.  The  overall  objective  for  the 
current  research  is  to  provide  an  improved  basis  for  making  damage  tolerance 
assessments  through  numerical  modeling  of  crack  tip  behavior,  including  the 
effects  of  plastic  or  other  residual  strains.  The  elastoplastic  BIE  method  is 
the  basis  for  the  current  effort. 

The  first  goal  of  the  originally  proposed  program  is  to  extend  an 
existing  planar  elastic  fracture  mechanics  analysis  based  on  the  BIE 
methodology  to  the  analysis  of  plastic  zones  around  cracks.  The  second 
proposed  goal  is  to  establish  fundamental  results  for  crack  tip  elastoplastic 
behavior,  based  on  a  numerical  and  analytical  study  of  the  elastoplastic  BIE 
formulation.  The  third  proposed  goal  is  to  establish  the  credibility  of  the 
elastoplastic  BIE  formulation  relative  to  the  finite  element  method  for 
refined  numerical  analysis  of  the  nonlinear  fracture  mechanics  problem,  and  to 
apply  the  capability  to  important  problems  of  fatigue  crack  growth  modeling 
for  advanced  aerospace  structures.  The  goal  for  the  second  year  of  the  effort 
is  to  extend  the  research  to  the  problem  of  modeling  crack  extension  under 
elastoplastic  conditions. 

The  specific  objectives  for  the  research  effort  as  stated  in  the  proposal 
are  as  follows.  Through  basic  research  in  the  elastoplastic  BIE  method, 
establish: 

1.  Cost  effective  and  highly  accurate  computational  method  for  planar 
fracture  mechanics  models  including  inelastic  effects  near  the  crack  tip. 

2.  Direct  comparisons  between  the  inelastic  BIE  fracture  mechanics  code 
and  an  advanced  finite  element  code. 

3.  Advanced  computational  procedure  for  investigating  history 
dependent  fatigue  cracK  growth  processes. 

JJ.  Feasibility  of  the  direct  numerical  analysis  of  elastic  and 
inelastic  response  to  growing  cracks:  the  results  of  this  capability  can  shed 
light  on  the  plastic  wake  phenomena  for  growing  fatigue  cracks. 
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All  first  year  goals  and  objectives  have  been  met.  New  understanding  of 
the  problems  associated  with  modeling  crack  growth  problems  has  been 
achieved.  The  next  section  summarizes  these  findings. 

RESEARCH  STATUS 
Review  of  the  Mathematics 

A  complete  treatment  of  the  elastic  formulation  for  the  Green's  function 
BIE  model  for  cracked  Dlanar  problems  is  given  by  Snyder  and  Cruse  [ 1 ]  and 
Cruse  [2].  The  full  development  of  the  elastopiastic  solution  is  given  by 
Cruse  and  Polch  [3].  The  following  summarizes  these  developments: 


ated  in  Figure 


The  basic  EIE  formulation  for  a  crack  problem,  as  illustr 
1 .  is  given  as  follows 


#  * 

C..u.(P)  +  Jt  (P,Q)u  (Q)ds  +  /T  (P,Q)u.(Q)ds 
J  s  J  J 

*  #  ® 

=  JU^P^t.^ds  +  JU^P.QJt.CQJds 


In  ( 1 )  the  u^,  tj_  terms  are  the  boundary  displacement  and  traction  vectors  for 
the  modeled  problem.  The  kernel  functions  (or  influence  functions)  U^*, 

T^j*,  are  mathematical  entities  giving  the  displacement  and  traction  that  are 
computed  on  S,  r  for  the  problem  of  an  infinite  body  loaded  at  p(x),  P(x)^  by 
a  set  of  unit  point  loads  in  each  coordinate  direction.  The  star  on  the 
kernel  functions  denotes  the  addition  to  the  point  load  solution  of  the  terms 
necessary  to  provide  for  a  traction  free  crack  at  a  specified  location  and 
orientation  in  the  geometry. 

The  use  of  a  Green's  function  for  spec  _il  geometries  is  well  developed  in 
potential  theory,  as  discussed  by  Greenberg  [4].  In  the  current  application 
we  seek  to  obtain  fracture  mechanics  solutions  for  the  case  of  traction  free 
cracks  in  finite  planar  bodies.  The  term  with  t^(Q)  for  Qer  in  (1)  is 
therefore  zero,  as  shown.  The  use  of  the  cracked  body  Green's  function 
results  in  the  traction  kernel  also  being  zero  on  the  crack,  viz. 

Tij*(P,Q)  =  0,Q  er,  also  as  shown  in  (1). 

Thus  (1)  constitutes  the  constraint  equation  that  must  be  satisfied  by 
u^,  tj^  on  the  uncracked  portion  of  the  surface.  This  equation  can  be  reduced 
to  solvable,  algebraic  form  through  the  use  of  suitable  approximations  to  the 


+Lower  case  p(x)is  an  interior  point:  upper  case  P(x)  is  a  boundary  point. 


boundary  data  uit  t^.  In  the  current  application  we  use  the  approximation  of 
piecewise  linear  interpolations  of  u,.,  ti  as  developed  by  Cruse  [2]. 

The  form  of  ( 1 )  for  the  interior  displacement  provides  a  means  of  direct 
computation  of  interior  strains,  stresses,  and  stress  intensity  factors. 

Simply  stated,  the  interior  quantities  depend  on  the  totality  of  boundary  data 
for  uit  t^  through  integration  of  these  quantities  together  with  appropriate 

kernel  functions  for  the  cracked  plane. 

Introduction  of  anelastic  strains  (e.g.,  residual  strains  due  to  welding, 

thermal  gradient  strains,  eiastopiastic  strains)  tilk  to  the  EIE  formulation  T* 
results  in  a  modification  to  (1) 


*(F,Q)ui(Q)cs 


;’J.i\p,Q)t.(Q)ds 

S'" 


-  Z*  (P,q)e^(q)dA 
<A>  J 


(2) 


The  addition  of  the  volumetric  (area  in  2D)  integral  in  (2)  is  seen  as  a 
correction  term  to  the  elastic  BIE,  (1).  The  kernel  function  z  . .  *  in  this 

*  i  ^ 

new  integral  consists  of  derivatives  of  the  elastic  displacement;  kernel  1*, 
and  its  form  differs  for  plane  stress  or  plane  strain. [5 j 


Equation  (2)  no  longer  provides  a  direct  means  for  computing  the  boundary 
data,  except  when  ei-H(q)  is  specified.  Thus,  for  eiastopiastic  response,  an 
additional  relationship  is  needed  to  compute  the  plastic  strains  for  (2).  The 
appropriate  equation  is  the  interior  strain  distribution,  as  written  by  Cruse 
and  Polch  [3]. 


"ij(p)  =  J*  Skij(p’Q)tk(Q)ds  +  JD*.  .(p,Q)uk(Q)ds 


S 


(3) 


This  equation  (3)  computes  the  interior  total  strain  in  terms  of  the 
boundary  data  and  the  interior  anelastic  strains.  For  eiastopiastic 

solutions,  the  unknown  data  ,  t^,  are  solved  for  incrementally  and 

equations  (2),  (3)  are  coupled  on  an  iterative  basis.  The  interior  anelastic 


strains  are  modeled  as  piecewise  constant  over  AA^area  segments  in  the  current 

study.  The  full  solution  algorithm  for  the  eiastopiastic  case  is  given  in 
Figure  2.  The  yield  criterion  has  to  be  satisfied,  giving  the  amount  of  total 
strain  that  is  plastic  at  each  load  level.  Tne  use  of  iteration  as  opposed  to 


The  dots  on  the  variables  denote  an  increment  in  the  variable. 


CONVERGED 


Figure  2.  Elastoplastic  Solution  Algorithm 


a  tangent  modulus  formulation  allows  us  to  precompute  all  of  the  elastic 
kernel  functions,  to  invert  one  of  these,  and  to  perform  all  of  the  ensuing 
numerics  as  matrix  multiplications. 


Stress  Intensity  Factor  Computations 

The  structure  of  (3)  has  been  investigated  by  Cruse  and  Polch  [3]  for 
interior  points  approaching  the  crack  tip.  It  was  found  that,  for  the 
elastoplastic  case  where  the  crack  tip  strains  can  exhibit  a  singularity  up  to 
1/o  (where  p  is  the  distance  from  the  crack  tip),  eq .  (3)  still  results  in 
convergent  integrals  ir.  (2),  (3).  However,  the  actua.  strength  of  the  plastic 
strain  singularity  is  a  function  of  the  work  hardening  (see  Hutchinson  [6]) 
and  car  only  be  inferred  from  the  resulting  strain  distr ioutions  after 
satisfying  tne  yield  criterion  implicit  in  (3). 


The  elastic  stress  intensity  factor  comoutation  for  the  EIE  formulation 


using  the  cracked  Green's  function  results  directly  from  the  elastic  version 
of  (3).  As  shown  by  Snyder  and  Cruse  [1],  the  kernels  in  the  two  boundary 
integrals  in  (3)  are  explicitly  dependent  on  the  inverse-square  root  of  the 
distance  of  p(x)  from  the  crack  tips  (+/-a).  Further,  for  the  nonsingular 
distribution  of  anelastic  strains  in  (3),  the  volumetric  kernel  has  the  same 
explicit  dependence.  Thus,  for  nonsingular,  anelastic  strains  we  obtain  the 
following  direct,  path  independent  evaluation  of  the  elastic  stress  intensity 
factors 


(K  ,K  )  =  -  JR  1  ’ 1 1  (Q)u  (Q)ds  «■  J L  I’II(Q)t  (Q)ds 

I  II  i  i  i  i 

s  s 

+  /M.  . 1 ’ 1 1 ( q ) e .  .A(q)dA 

1  j  1 J 

A 


The  first  two  terms  in  (4)  are  those  previously  used  by  Snyder  and  Cruse 
[1]  and  by  Stern,  et  al.  [7].  These  are  path  independent  integrals  which 
provide  a  simple  quadrature  for  computing  Kj,  Kjj  from  any  solution  for  u- ,  t^ 
on  a  path  around  the  crack,  but  excluding  the  crack. 

Equation  (4)  states  that  nonsingular,  anelastic  strains  modify  the 
elastic  K  j ,  Kjj  values  in  an  equally  simple  sense  of  quadrature  when  these 
quantities  are  specified  in  the  volume  (area).  Some  examples  of  this 
quadrature  for  a  notch  plasticity  problem  will  be  discussed  below. 

In  developing  eq.  (4),  it  was  assumed  that  the  anelastic  strains  were 
nonsingular,  thus  neglecting  the  crack  tip  elastoplastic  effects.  The 
additional  terms,  reflecting  the  higher  order  singular  behavior,  are 
represented  by  tne  incremental  elastoplastic  strain  portion  of  eq.  (3) 


(5) 
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As  discussed  by  Cruse  and  Polch  [3],  eq.  (5)  is  dimensionally  homogeneous  for 
any,  physical  singularity  in  plastic  strain  increment  as  p(x)+/-a,  but  the 
order  of  the  singularity  is  not  directly  soivacie  from  (5). 

A  considerable  number  of  technologically  important  problems  to  the 
aerospace  industry  are  associated  with  the  use  of  linear  elastic  fracture 
mechanics  parameters  (i.e.,  KT,  KT7;  for  problems  of  limited  or  localized 
plasticity.  These  include  predicting  KT  for  cracks  which  are  undergoing 
cyclic  plasticity  resulting  in  crack  closure  effects  on  spectrum  crack  growth, 
and  cracks  growing  ir  the  plastic  zone  of  a  bolthole  subject  to  high  loading 
or  prestressing. 

The  present  research  seeks  to  shed  light  on  some  of  these  problems  by 
presenting  a  stress  intensity  factor  computation  algorithm  that  can  directly 
and  unambiguously  model  these  kinds  of  limited  plasticity  effects.  For  such 
problems,  the  solution  given  in  (4)  is  to  be  used.  The  two  boundary  data 
integrals  in  (4)  reflect  the  plastic  strain  distribution  of  the  crack  tip,  as 
well  as  other  ;nelastic  strains  through  the  volume  integral  in  (2).  Secondly, 
the  prior  plasticity  of  a  notch  will  affect  Kj,  Kjj  through  the.  volume 
integral  of  those  nonsingular  strains  in  (4).  The  resulting  values  of  Kj,  Kjj 
are  the  plasticity  corrected  elastic  stress  intensity  factors  which  define  the 
strength  of  the  elastic  singularity  which  dominates  the  plastic  singularity. 
The  use  of  this  approach  is  obviously  limited  to  crack  tip  plastic  zones  which 
are  contained  within  the  field  of  the  elastic  singularity. 

Numerical  Solution  Algorithm 


Application  of  the  appropriate  interpolations  to  the  data  in  eqs.  (2)  and 
(3)  reduces  the  integrals  to  algebraic  form.  In  general,  the  boundary 
solution  involves  an  equal  number  of  known  (applied)  boundary  data  and  unknown 
data.  Letting  the  unknown  data  be  given  by  {*},  the  product  of  the  known  data 
and  its  coefficient  matrix  terms  by  {^},  and  the  coefficient  of  the  piecewise 
constant  plastic  strains  by  [E],  we  obtain  from  (2) 

[A] {x}  =  {y}  +  [E]{eP}  (6) 

Similarly,  taking  [S]  and  [D]  to  be  the  elastic  coefficient  arrays  of 
the  boundary  data,  and  [G]  to  be  the  elastic  coefficient  array  for  the 
plastic  strain,  then  eq.  (3)  becomes 

{ eT }  =  [S]{t}  -  [D] {u}  -  [G] {ep}  (7) 

The  strain  superscripts  in  (6)  and  (7)  refer  to  total  (elastic  plus  plastic) 
and  plastic  values,  while  the  dots  imply  that  all  of  the  variables  are  to  be 
interpreted  in  terms  of  their  incremental  evaluation. 
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The  present  BEM  algorithm  makes  use  of  the  Huber-Mises-Hencky  yield 
condition  and  associated  flow  rule.  Elastic,  perfectly  plastic  material 
response  was  modeled  throughout  this  study,  but  the  code  allows  for  a  multi- 
piecewise-linear  definition  of  a  general  stress-strain  curve. 

Following  the  approach  adopted  in  the  ADINA  code,  each  increment  in  total 
strain  is  divided  into  subincrements  (Bathe  [8]).  The  number  of  subincrements 
is  selected  to  minimize  the  error  in  the  deviatoric  stress  change  within  the 
strain  increment.  The  stress  increment  for  a  given  iterate  of  increment  in 
plastic  strain  is  then  obtained  by  an  Euler  forward  integration  of  the  flow 
rule,  according  to 

T 

CURRENT 

LOAD 


{c  }  =  {o  ’  *  [C  ]d  * 

CURRENT  PREV  J  (8) 

LOAD  LOAD  T 

PREV 

LOAD 

In  eq.-  (8),  the  elastoplastic  matrix  relating  the  subincrements  in  stress 
and  total  strain  is  given  for  plane  strain  by 
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The  current  state  of  deviatoric  stress,  S^,  and  second  invariant  of 
deviatoric  stress,  J2,  is  updated  within  the  subincremental  integration  of 
(8).  The  tangent  modulus,  H,  is  taken  as  the  slope  of  the  effective  stress- 
effective  plastic  strain  curve,  at  the  current  level  of  effective  stress. 

Figure  2  summarizes  the  current  iteration  algorithm  for  the  solution  of 
eqs.  (6,7).  The  coefficient  arrays  [A],  [S],  and  [D]  depend  solely  on  the 
elastic  constants  of  the  material  and  the  boundary  shape.  Thus,  they  are 
computed  once  and  stored.  The  [A]  matrix  is  inverted  prior  to  storage.  The 
interior  arrays  [E]  and  [G]  are  also  dependent  solely  on  the  elastic  constants 
and  the  interior  element  modeling.  These  are  also  computed  once  and  stored. 
Note  again  that  only  that  portion  of  the  interior  expected  to  be  inelastic 
need  be  modeled.  The  expense  of  generating  [E]  and  [G]  for  crack  problems 
dictates  that  such  limited  volumetric  modeling  be  employed. 

In  the  first  iteration  at  a  given  load  step,  the  plastic  strain  increment 
in  Figure  2  is  taken  from  the  last  load  step.  The  boundary  solution  then 
responds,  in  an  elastic  manner,  to  the  increase  in  loading.  Estimated 
interior  total  strains  are  then  calculated.  Based  on  the  new  total  strain 
increment,  the  interior  stresses  and  plastic  strains  are  computed  based  on 
satisfying  the  yield  condition  through  eq.  (8).  The  plastic  strain  increment 
is  then  updated  in  both  eqs.  (6,7)  for  a  recalculation  of  the  boundary  and 
interior  solutions. 


Absolute  convergence  of  the  strain  solution  within  each  element  is 
required  for  the  iteration  process  usee.  That  is,  the  maximum  difference 


between  successive  iterates  of  the  plastic  strain  correction  term  (second  term 
on  the  right  hand  side  of  eq.  (6))  is  not  allowed  to  exceed  a  user-specified 
tolerance.  This  tolerance  has  been  selected  on  the  basis  of  its  ability  to 
relate  directly  to  the  amount  of  the  displacement  increment.  A  number  of 
numerical  experiments  with  tolerances  ranging  over  10“^  to  10*9  were  conducted 
to  test  the  sensitivity  of  the  results  to  this  value.  It  was  found  that  the 
errors  in  the  displacements,  for  a  simple  uniform  stress  test  case,  were  of 
the  order  of  the  tolerances  specified.  A  decrease  in  the  tolerance  by  an 
order  of  magnitude  generally  resulted  in  a  doubling  of  the  number  of 
iterations  required  to  achieve  convergence.  A  value  of  10“^  was  used  for  the 
notch  problem  and  a  value  of  10"^  was  used  for  the  fracture  mechanics  problem, 
in  order  to  account  for  the  higher  strain  gradient. 

Numerical  Results 

The  computer  program  has  been  verified  on  two  example  problems.  The 
first  is  a  plate  with  perforations,  previously  solved  by  Haward  and  Owen  [9] 
using  finite  elements  and  resolved  by  Telles  [10]  using  the  BEM.  This  example 
served  the  basic  purpose  of  validating  the  current  code  and  provided  a  basis 
for  some  numerical  experimentation.  The  second  problem  is  a  fracture 
mechanics  problem  of  a  center  cracked  plate  loaded  in  tension.  The  plastic 
strain  results  are  compared  to  ADINA  results  using  a  singular  finite  element 
model . 

The  geometry  for  the  first  problem  is  shown  in  Figure  3.  Plane  strain 
conditions  are  applied  for  all  three  of  the  analyses  and  the  material  is  taken 
to  be  elastic-perfectly  plastic.  The  appropriate  constants  are  E  =  42.  x 
10^MN/m2,  o  =  105.  MN/m^,  v  =  0.33.  The  one  loading  condition  considered 
was  uniaxial  tension,  applied  by  prescribing  displacements  at  the  edges  of  the 
plate  section.  The  piecewise  linear  plastic  strain  BEM  mesh  of  Telles  is 
shown  in  Figure  3;  the  FEM  quadratic  isoparametric  element  mesh  used  by  Haward 
and  Owen  is  shown  in  Figure  4.  The  current  BEM  mesh,  using  piecewise  constant 
plastic  strains,  is  shown  in  Figure  5. 

Figure  6  plots  the  numerical  results  in  terms  of  the  amount  of  force 
required  versus  the  applied  displacements.  Table  1  summarizes  the  numerical 
force-displacement  data.  All  three  model  results  show  excellent  agreement, 
given  the  disparity  in  modeling  strategies.  The  predicted  limit  load  for  the 
current  study  differs  from  the  other  two  by  less  than  2%.  The  difference  is 
attributed  to  the  use  of  constant  strain  elements.  Limit  load  is  obtained 
when  the  centroidal  value  of  stress  in  the  last  ligament  element  yields,  a 
condition  that  will  occur  below  the  load  for  yielding  the  last  physical 
ligament  ahead  of  the  notch. 

The  current  BEM  code  was  tested  for  a  range  of  load  increments  in  a 
deliberate  attempt  to  create  numerical  instability.  The  numerical  results 
plotted  in  Figure  7  were  generated  using  an  incrementation  scheme  resulting  in 
a  single  element  yielding  at  a  time.  The  BEM  results  required  about  twenty 
iterations  per  load  step  to  fully  converge.  The  worst  case  was  one  load  step 
to  the  maximum  displacement.  The  solution  converged  in  45  iterations  and 
agreed  with  the  other  limit  load  results  within  0.2 %.  The  maximum  deviation 
in  calculated  plastic  strains  was  10K  in  the  last  element  to  yield. 


Figure  3.  Perforated  Polystyrene  Plate  and  BIE  Discretization 

Used  by  Telles  (1983) 
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Table  1. 

Numerical  Results  of  the  Polystyrene  Plate  Problem 


Load  Case 


PCP  -  Uniaxial  Stretch 
Load  %  Displacement  Force 


■  (x  10-2)  F/(a+d) 


10 

.3576 

11 

.3259 

12 

.2459 

13 

.1047 

13 

.8412 

14 

.4610 

14 

.7765 

14 

.9091 

15 

.0009 

15 

.0650 

15 

.1112 

15. 

1412 

15. 

1412 

40.820 


50.003 

50.217 


The  ultimate  strength  =  50.471  MN/m 
a,  d  defined  in  Figure  2 


The  second  example  is  a  center-cracked  plate  loaded  in  tension.  The 
total  width  of  the  plate  is  8  units,  with  a  crack  size  of  2  units.  One 
quarter  of  the  geometry  was  modeled  using  ADINA,  as  shown  in  Figures  8  and 
9.  Extremely  fine  resolution  of  the  crack  tip  elements  was  taken  in  order  to 
minimize  the  error  in  the  finite  element  solution.  The  elastic  stress 
intensity  factor  for  this  finite  element  model,  using  the  crack  opening 
displacement  at  the  quarter-point  node,  was  in  error  relative  to  handbook 
values  by  about  2%. 

The  BEM  mesh  corresponding  to  the  local  finite  element  scale  is  shown  in 
Figure  10.  The  elastic  3EK  stress  intensity  factor  results  were 
indistinguishable  from  the  handbook  results.  The  FEM/BEM  meshes  were  selected 
so  as  to  provide  about  three  decades  of  plotting  data  in  terms  of  crack  tip 
distance.  The  maximum  size  of  the  plastic  zone  was  limited  solely  for 
convenience  in  the  current  study. 

Plane  strain  conditions  were  used _ for  both  of  the  models.  The  elastoplast 
material  constant  used  were  E=2.037- 103MN/m2 ,  0^=3. ^52* 102MN/m: ,  and  v=0.27. 

The  ADINA  crack  tip  model  used  quadratic,  isoparametric  finite  elements 
with  nine  interior  strain  integration  points.  The  BEM  model  used  constant 
strain  triangles  throughout.  As  noted  above,  both  meshes  were  identical  in 
the  crack  tip  region.  The  ADINA  model  used  one  layer  of  collapsed  quadratic 
elements  adjacent  to  the  crack  tip.  This  approach  induces  a  (1/r)  type  of 
singularity  in  the  displacement  gradient  within  this  first  layer  of 
elements.  No  singularity  modeling  is  used  in  the  BEM  plastic  strain 
distribution. 

Loading  history  was  identical  for  both  models  and  spans  the  range  of  load 
factors  of  0.0310  to  0.2075.  A  value  of  1.0  corresponds  to  yielding  of  the 
whole  plate.  A  total  of  68  load  steps  was  used  for  both  models.  The  load 
steps  satisfy  the  conditions  of  Larsson  and  Carlsson  [11].  Simply  stated, 
these  conditions  require  that  at  most  one  element  becomes  plastic  at  each  load 
increment,  and  that  the  load  increment  should  be  smaller  than  1%  of  the  load 
corresponding  to  KImax=°y The  range  of  load  factors  has  been  chosen  as 

the  range  to  go  from  yielding  the  innermost  element  to  yielding  the  outermost 
element . 

ADINA  failed  to  converge  for  the  first  step  until  the  stiffness 
reformulation  (BFGS)  procedure  was  used.  After  the  first  load  step  (requiring 
20  iterations)  the  ADINA  algorithm  with  reformulation  generally  converged  with 
five  iterations.  The  BEM  algorithm,  using  elastic  "stiffnesses,"  converged  in 
ten  to  fifty  increments  at  each  load  step  with  the  higher  numbers  occurring  at 
the  higher  load  levels.  The  total  computer  time  for  the  two  models  was 
essentially  the  same,  although  the  BEM  calculations  are  cheaper  per  load 
step.  A  higher  final  load  level  or  cyclic  loading  would  yield  a  benefit  to 
the  BEM  model,  even  though  the  current  BEM  code  is  not  yet  optimized  for  these 
calculations. 

The  crack  tip  plastic  strain  distribution  results  are  shown  in  Figure  11 
for  two  of  the  computed  load  levels.  The  data  are  taken  from  points 
distributed  near,  but  not  on,  a  line  at  an  angle  of  about  85°  to  the  plane  of 
the  crack.  This  angle  corresponds  to  the  line  of  maximum  equivalent  elastic 
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strain.  The  jaggedness  of  the  curves  is  mostly  due  to  the  use  of  triangular 
elements,  as  well  as  to  the  points  having  different  angular  locations.  The 
data  is  plotted  in  terms  of  the  centroidal  value  of  plastic  strain.  The 
innermost  row  of  finite  elements  has  three  sampling  points  radially,  which 
accounts  for  the  smaller  radius  plotted  for  these  results.  The  numerical 
results  from  ADINA  show  a  tendency  for  a  spurious  peak  in  plastic  strain 
increments  in  the  second  row  of  elements.  This  peak  is  no  doubt  induced  by 
the  lack  of  a  singularity-transition  element  in  the  current  study. 

It  is  significant  to  find  that  the  numerical  results  are  in  such  good 
agreement.  This  confirms  the  accuracy  of  the  BEM  algorithm  for  piecewise 
constant  plastic  strains.  The  BEM  results  do  not  show  the  strength  of  the 
plastic  singularity  as  strongly  in  the  first  row  of  elements  as  do  the  finite 
element  results,  with  the  imbedded  1/r  singularity  in  displacement  gradient. 
However,  both  sets  of  results  strongly  indicate  that  the  plastic  strain  for 
localized  plasticity  possesses  the  same  1/r  singularity  that  is  associated 
with  fully  developed  plasticity  for  the  case  of  zero  strain  hardening. 

Clearly,  the  presence  of  the  underlying  elastic  singularity  field  plays  an 
important  role  in  enhancing  the  modeling  accuracy  for  crack  tip  plasticity. 

Figures  12  and  13  show  the  progressive  development  of  the  plastic  zone  up 
to  the  maximum  modeled  load.  It  is  to  be  emphasized  that  the  current  study 
was  intended  to  confirm  the  accuracy  of  the  new  BEM  algorithm  for 
elastoplastic  fracture  mechanics  analysis,  and  not  study  extensive  plasticity 
response  at  the  crack  tip.  For  this  reason  the  current  results  were  not 
carried  beyond  the  load  level  shown.  There  is  no  inherent  limit  to  the  load 
level  that  can  be  modeled  with  this  BEM  algorithm. 

The  next  problem  was  selected  to  validate  the  stress  intensity  factor 
algorithm  for  prior  plasticity  for  any  residual  or  thermal  strain  field.  The 
geometry  selected  is  a  simple  tension  specimen  with  the  boundary  and  internal 
mesh  shown  in  Figure  14.  The  mesh  arrangement  was  selected  solely  for 
convenience,  as  it  is  used  as  a  portion  of  a  later  mesh.  The  specimen  was 
loaded  to  110%  of  the  yield  stress  for  a  bilinear  material  response.  This 
induced  a  uniform  plastic  strain  throughout  the  specimen. 

The  next  step  in  the  validation  of  eq.  (4)  was  to  introduce  a  crack,  done 
along  the  bottom  of  the  mesh  as  shown.  The  residual  boundary  solution 
corresponding  to  the  residual  internal  strains  is  computed  for  the  cracked 
case  by  eq.  (6).  Next,  the  internal  strains  for  the  residual  boundary  and 
internal  variables  are  computed  from  eq.  (7).  In  the  case  of  the  test 
problem,  eq.  (6)  produced  the  uniform  displacements  compatible  with  uniform 
residual  strains,  eq.  (7)  computed  internal  strains  equal  to  the  residual 
strains. 

The  elastic  stress  intensity  factor  for  the  problem  was  then  computed  for 
the  residual  boundary  terms  computed  from  (6).  If  there  were  further  changes 
in  the  residual  strains  due  to  unloading  plasticity,  these  would  modify  the 
elastic  strain  intensity  factor  through  the  appropriate  term  in  eq.  (4).  As 
required  for  this  simple  case,  the  residual  stress  intensity  factor  was 
zero.  The  residual  values  in  the  first  and  third  terms  in  eq.  (4)  cancel  each 
other  to  within  computer  accuracy. 
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Figure  15  shows  the  final  application  problem  for  this  study.  The  model 
includes  a  simplified  internal  mesh  for  modeling  plasticity  at  a  hole  in 
tension.  The  plate  is  two  units  long  by  one  unit  wide  (W).  The  hole  radius 
(R)  is  0.25W.  Local  meshes  were  used  to  model  crack  tips  at  0.05R  and  at 
0.5R.  The  elastic  KT  for  this  mesh  is  computed  to  be  3.29  as  compared  to  3.24 
from  Peterson  [12].  The  elastic  stress  intensity  factors  for  the  two  crack 
lengths  are  13 %  higher  than  the  values  given  by  Rooke  and  Cartwright  [13], 
which  is  due  to  the  finite  width  of  the  specimen. 

Figure  16  shows  the  progressive  generation  of  the  plastic  zone  at  the 
hole.  The  material  is  modeled  in  plane  strain  with  a  yield  stress  (perfect 
plasticity)  of  50  ksi.  The  plate  is  loaded  to  a  level  of  44.6  ksi  applied  at 
the  end  of  the  plate.  This  level  of  loading  was  exaggerated  in  order  to 
generate  a  plastic  zone  that  encompassed  both  crack  lengths.  In  fact,  the 
notch  was  found  to  undergo  substantial  reverse  yielding  at  unloading.  After 
the  plastic  strains  are  computed  for  the  notch  alone,  the  crack  is  introduced 
through  eqs.  (6,7).  The  boundary  solution  computes  crack  face  closure 
corresponding  to  the  residual  strains.  With  application  of  the  residual 
strain  and  boundary  terms,  the  internal  strain  algorithm  predicts  localized 
reversed  plastic  flow  at  the  crack  tip.  Elastic  stress  intensity  factors  are 
computed  for  the  residual  terms,  without  the  unloading  corresponding  to 
introduction  of  the  crack,  and  for  the  residual  terms  with  full  unloading  of 
the  crack  surface. 

Table  2  summarizes  the  numerical  results.  The  engineering  approach  often 
used  for  cracks  at  notches  with  prior  plasticity  is  to  take  the  residual 
strains  (and  resulting  stresses)  prior  to  unloading  and  to  use  these  as  pseudo 
tractions  on  the  crack  faces.  This  approach  corresponds  to  the  second  set  of 
results  in  Table  2,  for  zero  crack  surface  unloading.  The  results  show  that 
the  residual  strains  produce  crack  closure  (negative  stress  intensity  factor) 
and  a  corresponding  reduction  in  the  maximum  stress  intensity  factor  at 
load.  The  effect  is,  of  course,  more  pronounced  for  the  small  flaw  as  it  is 
more  completely  buried  by  the  prior  plastic  zone. 


Table  2. 

Stress  Intensity  Factor  Results  (KI/o/ttS) 


a/P  =  C.C5 

a'R  -  0.5 

Elastic 

-  “lax.  load 

-  Zero  load 

0 

2.070 

0 

Plastic  (Zero  Crack 

Surface  Unloading) 

-  Max.  load 

-  Zero  load 

1  .«83a 
-1.923s 

1 ,886a 
-0. 18*.b 

Piastic  (Full  Crack 

Surface  Unloading) 

-  Max ,  load 

-  Zero  load 

O.B«?a 
-2 . 568s 

»  ?37® 
-0. ;82s 

Note  a:  Maximum  load  values 
Elastic  *  2ero  load 

(elastically)  = 
values 

Note  b:  Negative  values  of  Kj  imply  crac*  Closure 
at  positive  load 


The  more  correct  result  for  the  stress  intensity  factor  calculations  is 
provided  by  the  third  set  of  results  in  Table  2  with  full  crack  surface 
unloading.  The  point  to  be  made  is  the  significant  difference  between  the 
results  assuming  no  crack  face  unloading  of  the  residual  stresses,  and  those 
for  full  unloading.  The  usual  engineering  approaches  to  modeling  crack  growth 
in  a  notch  with  plasticity  [14,  15]  do  not  account  for  this  unloading,  which 
occurs  at  positive  applied  loads  (crack  open).  The  effect  is  most  significant 
for  short  cracks,  which  generally  consume  most  of  the  crack  growth  life. 

Tne  same  discrepancy  between  the  engineering  and  more  correct  results, 
including  unloading  effects,  will  occur  for  problems  of  residual  welding 
strains  and  thermal  gradient  strains.  The  new  algorithm  presented  herein  is 
offered  as  an  effective  means  for  correcting  the  engineering  approach  when  the 
plastic  zone  can  be  estimated  or  computed,  if  not  by  BIE,  then  by  finite 
element  procedures. 

RESEARCH  COMMUNICATION 

The  research  to  date  has  involved  three  individuals.  The  program  manager 
is  Dr.  T.  A.  Cruse.  Dr.  Arje  Nachman  was  a  co-principal  investigator  on  the 
effort  for  a  few  months  prior  to  his  taking  a  position  with  the  Hampton 
Institute,  Hampton,  Virginia.  Dr.  E.  2.  Polch,  who  was  developing  the 
computer  algorithms  for  the  current  research,  has  been  serving  as  co-principal 
investigator  since  chat  time. 

The  research  results  to  date  have  been  incorporated  in  several 
manuscripts.  These  include  a  chapter  in  a  major  reference  book  on  BIE,  two 
journal  submittals,  and  two  conference  proceedings.  The  following  list 
summarizes  these  submittals: 

1.  "Fracture  Mechanics,"  T.  A.  Cruse,  Boundary  Element  Methods  in 
Mechanics,  book  in  series  Computational  Methods  in  Mechanics,  edited  by  D.  E. 
Beskos,  Elsevier  Science  Publishers  B.V.,  Amsterdam,  to  be  published. 

2.  "Elastoplastic  BIE  Analysis  of  Cracked  Plates  and  Related  Problems, 
Part  1:  Formulation,"  T.  A.  Cruse  and  E.  Z.  Polch,  submitted  to  International 
Journal  for  Numerical  Methods  in  Eneineerine. 


3.  "Elastoplastic  BIE  Analysis  of  Cracked  Plates  and  Related  Problems, 
Part  2:  Numerical  Results,"  T.  A.  Cruse  and  E.  Z.  Polch,  submitted  to 
International  Journal  for  Numerical  Methods  in  Engineering. 

4.  "Advanced  Algorithms  for  Fracture  Mechanics  Analysis  in  Two  and  Three 
Dimensions,"  T.  A.  Cruse  and  E.  Z.  Polch,  2nd  International  Conference  on 
Variational  Methods  in  Engineering,  Southampton,  England,  July  17-19,  1985. 

5.  "BIE  Analysis  of  Crack  Tip  Plastic  Zones,"  T.  A.  Cruse  and  E.  Z. 
Polch,  AIAA  26th  Structures,  Structural  Dynamics,  and  Materials  Conference, 
Orlando,  Florida,  April  15-17,  1985. 

The  paper  from  the  26th  AIAA  Structures,  Dynamics,  and  Materials 
conference  is  currently  under  revision  to  add  some  additional  information 
relative  to  the  calculation  of  crack  tip  stress  intensity  factors,  including 
the  effects  of  plasticity.  The  paper  will  be  submitted  shortly  to  the 


International  Journal  of  Fracture  and  will  also  be  co-authored  by  Cruse  and 
Polch.  The  efforts  from  the  second  year  of  the  contract  will  also  result  in 
several  publication  submittals.  At  this  time  it  is  probable  that  these 
articles  will  be  submitted  for  publication  in  the  Engineering  Fracture 
Mechanics  Journal  and  the  International  Journal  of  Fracture.  Some  specific 
aspects  of  the  current  research  will  be  included  in  further  review  type 
articles  by  Cruse,  and  in  various  conference  proceedings. 

Dr.  Cruse  has  made  two  trips  associated  with  other  business  activities 
which  have  resulted  in  specific  consultation  on  the  problems  of  crack  growth 
modeling  with  elastoplastic  models.  The  current  research  effort  was  aiscussec 
in  detail  with  Mr.  James  Rudd  of  the  Air  Force  Flight  Dynamics  Laboratory  on 
and  with  Dr.  James  Newman  of  NASA  (Langley  Research  Center).  Botn  individuals 
expressed  interest  in  this  work  and  will  be  kept  informed  of  the  research 
progress.  A  brief  telecon  was  held  with  Dr.  Ted  Nicholas  to  inform  him  of  the 
principal  results  of  the  work  and  to  inform  him  of  the  potential  for 
application  of  the  research  to  the  small  flaw  problem. 
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